Quantum quenches and off-equilibrium dynamical transition in the 
infinite-dimensional Bose-Hubbard model 
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We study the ofT-equilibrium dynamics of the infinite dimensional Bose Hubbard Model after a 
quantum quench. The dynamics can be analyzed exactly by mapping it to an effective Newtonian 
evolution. For integer filling, we find a dynamical transition separating regimes of small and large 
quantum quenches starting from the superfiuid state. This transition is very similar to the one 
found for the fermionic Hubbard model by mean field approximations. 



Significant advances in the field of ultra cold atoms 
have allowed one to engineer quantum many-body sys- 
tems in almost perfect isolation from the environment. 
Thanks to the abihty to rapidly tunc different param- 
eters, e.g. the interaction strength between the atoms 
or the creation of controlled excitations, the realm of 
non-equilibrium many body physics of (almost) isolated 
quantum systems has thus been accessed and can now be 
studied experimentally. For example, Greiner et al. [1| 
studied the dynamics of interacting bosons loaded on an 
optical lattice. The physics of this system is well cap- 
tured by the Bose-Hubbard model. By changing the in- 
tensity of the lasers one can effectively tune the param- 
eters in the corresponding Bose-Hubbard model. Rapid 
changes induce interesting non-equilibrium dynamics [l| . 
The activity in this field is booming: several new ex- 
periments have been performed, including on fermionic 
systems l2,j3| ; questions about thermalization [J, Isf , its 
absence [6|-|8[ , quantum dynamical phase transitions out 
of equilibrium [9|, llO| are currently addressed. 
A protocol inducing an off-equilibrium dynamics, which 
has received a lot of attention recently, is the so called 
quantum quench. It corresponds to preparing the system 
in the ground state of the Hamiltonian Hi, to changing 
suddenly at time t = a parameter of the Hamiltonian, 
for example the interaction strength, and then letting 
the system evolve with the new Hamiltonian Hf. Sev- 
eral studies have been performed for the Bose Hubbard 
model, which as discussed above is relevant for experi- 
ments. There have been numerical analysis of one di- 
mensional systems by exact diagonalization and t-DMRG 
[J, a LZI, ISl- Saddle point approximation [ll[, Gross- 
Pitaevskii equations [12| and Gutzwiller approximation 
[ij, [l^ have been used to analyze higher dimensional 
and realistic cases. The fermionic Hubbard model has 
been also studied by mean field theories recently [9|, llO| . 
In this work we present a complete analysis of quantum 
quenches in the Bose Hubbard model (BHM) in the limit 
of infinite dimensions. The advantage of this limit is that 
the model can then be analyzed exactly even out of equi- 
librium. Its solution at equilibrium played an important 
role in determining the phase diagram and the properties 
of the Mott-supcrfluid quantum phase transition of the 



three dimensional BHM [l3[ . Studying its off-equilibrium 
dynamics is therefore a natural route to follow. We will 
discuss in the conclusion the limitations of this approach 
and possible extensions. To obtain a well defined infinite 
dimensional limit one has to scale the hopping amplitude 
as one over the dimension d[l6|. A complementary but in 
the bosonic case identical procedure [? ], which we will 
follow for simplicity, consists in focusing from the start 
on the BHM defined on a completely connected lattice. 
The corresponding Hamiltonian reads: 
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where 6j , bi are the bosonic creation and annihilation 
operators, fii = blbi the occupation operator and V the 
total number of sites. In the following we take J = 1 and 
measure U in units of J and the time in units of J/fi. We 
shall study off-equilibrium dynamics induced by quan- 
tum quenches corresponding to a sudden change of the 
interaction strength from Ui to Uf at t = 0. Since H is 
invariant under any permutation of sites, all eigenstatcs 
can be classified in terms of the corresponding symmetry 
classes. In particular the ground state, whether Mott or 
superfiuid, corresponds to a completely site permutation 
symmetric wavefunction. Since also the time-dependent 
wavefunction remains completely symmetric after the 
quench, one can restrict the analysis to the subspace of 
completely symmetric states. It is easy to convince one- 
self that these states can be parametrized by the frac- 
tion xo, xi, X2, ... of sites with 0, 1, 2, . . . bosons and that 
they correspond to the flat normalized sum of all Fock 
states characterized by Vxi sites with i bosons per site. 
In order to simplify the presentation, let us first focus 
on the simplified model where maximum two bosons per 
site are allowed, ni, = 2. We shall discuss later the gen- 
eralization to any value of n?,. Since xq -|- xi -I- a;2 = 1 for 
rib = 2 and because the number of particles T^(a;i -1-2x2) is 
conserved by the dynamics, a generic symmetric state is 
identified by xi only, and can be denoted \xi) (henceforth 
wc will drop the subindex 1). The evolution of the wave- 
function 1-0) = '^^'>Px\x) is determined by the equation 

{x\idtJ2x' ^x'\^') = {x\HJ2x''^^'\^')- ^^ this model, all 
matrix elements {x\H\x') are zero except {x\H\x ± 2/V) 
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FIG. 1. Graphical solution for the value of p at the turning 
points. The trajectories are full lines, and the position at 
t = is indicated by a circle for the three trajectories A, B 
and C. In case A it is impossible to have a turning point at 
p — -K 12. Case B corresponds to the dynamical transition, 
and C to unbounded evolution of p. A, B and G are plotted 
in Fig. [3 



and the diagonal term {x\il\x)\ the former corresponds 
to the physical process of one boson jumping from one 
site to another. The resulting Schrodinger equation for 
■\l)x reads : 



idtiix = D{x)i)j; - W{x) i^tpx+2/V + '(Px^2/v) 

D{x)- 2W{x)cosh(2dx/V)\ijJx (2) 

D{x)- 2W{x)cos{2p))'i{jx 



where W{x) = x[{2-x~n){n-x)/2]^^'^, D{x) = U(n- 
x)/2— x{2 + n — 'ix)/2, n is the number of bosons per site 
and subleading contributions in 1/V have been dropped. 
The initial wavefunction, which is the ground state at 
couphng Ui, is a wave packet of width I/a/F, see [l8| and 
below. Since 1/V plays the role of 7i in ([2]) the thermody- 
namic limit corresponds to the classical limit. As a con- 
sequence, the time evolution of the average particle posi- 
tion x{t) — (x) and momentum p(t) = (p) = {~idx/V), 
is given by the Newton equations for the Hamiltonian 
H = D[x) — 2W{x) cos(2p), where x{t) and p{t) are clas- 
sical canonical variables. The validity of this argument 
can be thoroughly established by a direct analysis jlSf . 
In particular, one can show that on timescales less than 
y/V, ipxit) ^ e^p[V{x - x{t)f/2a{tf + iVp{t)x\, i.e. it 
is a sharp wave-packet, centered at x{t), of width of or- 
der l/yfV and has a very fast oscillating phase e^^pW^^ 
where x(t) and p{t) are the classical canonical variables 
defined above. In the following we will repeatedly make 
use of this mapping to a classical swtem. Similar map- 
pings have been recently used in 10|, ll9|, l20[ . The first 
useful consequence is that the ground state is obtained 
minimizing H with respect to p and x ; the corresponding 
p is actually always zero; in consequence the ground state 
is obtained by the value of x minimizing D[x) — 2W{x). 
The phase diagram is similar to the one derived by Fisher 
et al. [16| except that there is only one Mott lobe cor- 
responding to n = 1. As we shall see, it is at integer 
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FIG. 2. Evolution as a function of time of x, Ivl'oP (top panels) 
and p (bottom panels) increasing the amplitude of the quench. 
Uf — 3.33 is kept fixed and several Ui are considered, with 
A, B, G corresponding to Ui = {1.62,0.838,0.1} (unlike in 
the text where Uf is varied). The case B correspond to the 
transition point. 



filling, n ~ 1, where the Mott state exists, that the 
off-equilibrium dynamics is most interesting. We shall 
consider this case first, for which the ground state (GS) 
corresponds to 
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with f/c = 3 + 2V2. In this simple model, the con- 
densate fraction |*oP is simply equal to y^^^^jbjbi, 
which up to a sign coincides with the average value of 
the (intensive) kinetic energy. This can be easily ob- 
tained subtracting the average value of the interaction 
term to the total energy and reads, for the ground state, 
|*o|' = a:Gs(l-a^Gs)C/c/2. 

Let us now consider quenches starting from a superfluid 
ground state and increasing the value of U from Ui to 
Uf . A small increase of U leads to oscillations of x and p 
as can be verified analytically and checked numerically, 
see Fig [2j\. The turning points of x{t), determined by 
X — 4:W{x) sm{2p) = 0, correspond to p = 0. Actu- 
ally there would be the possibility to have p ~ mr/2 
too. However, a p starting from zero and reaching the 
value 7r/2 would imply, by energy conservation, a value 
of X at the turning point such that E = D{x) + 2W{x), 
where E is the energy after the quench. This equation 
has no solution for small quenches as shown graphically 
in Fig. [H see case A. It starts to have a solution for 
larger quenches, when E becomes positive (cases B and 
C in Fig. [1]). Actually Ed = corresponds to a dy- 
namical transition: for E < Ed the momentum p{t) is 
bounded, whereas for E > Ed it grows infinitely large. 
The condition £' = depends on Ui,Uf. One finds 
that for a given Ui, the corresponding critical value is 
Uf{Ui) = {Ui + Uc)/2. Approaching Uf the period of os- 



cillation increases and diverges as r 
where c 



''H\Uf-uf\), 

y/{Uc~Uf)iUf-l/Uc). Fig. M shows the 
typical time-evolution of cc, |^oP andp for the three cases 
A, B and C. Exactly at t/J? the system relaxes exponen- 
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FIG. 3. a) Evolution of time average (continuous line) and 
microcanonical average (dashed line) of (|*o|^) as a function 
of Uf for Ui = 0. b) Dynamical phase diagram for the model 
with maximum two bosons per site. Quenches from the Mott 
phase are not considered. Quenches from the superfluid phase 
are oscillating and similar to A or C. The dynamical transition 
separating the two is displayed as a dashed line, it meets the 
Mott phase at Uf = Uc- The phase diagram for the case of 
more than two boson per site is qualitatively similar. 
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FIG. 4. a) As Fig. [2]but for n^ = 3 and Ui = 1. Left and 
right: Uf = 2.5 and 3.29, respectively below IE < 0) and 
above {E > 0) the dynamical transition at Uf = 3.21. b) 
Variation of (I'I'oP) as a function of Uf for n = 1, [/^ = 3, 
with Ub = 2, 3, 4 and 5. The plot of ni, = 2 is shifted of 0.025 
along the Uf axis for comparison. 



tially to the Mott state with a rate c. Approaching the 
transition the system spends most of the time close to the 
Mott state and therefore the time averaged condensate 
fraction vanishes at Vi in a singular way, proportional to 
1/r. This singularity is related to the fact that the Mott 
state is 'absorbing': classical trajectories falling into it 
cannot escape, and the period t diverges when approach- 
ing ui. Conversely, trajectories starting from the Mott 
state remain stuck to .t = 1 on large times, t ~ log V . 
This is, however, a peculiarity of the infinite dimensional 
limit; for a finite dimensional system, spatial fluctuations 
will drive the system away from the Mott state [25|, |26[ . 
In Fig. 13^, as an example of singular behavior, we show 
l^oP as a function of Uf for quenches starting from the 
non interacting case C/^ = 0. Moreover, we compare 
IvpQp to its microcanonical average at the same energy. 
Clearly, the system is not thermalized. At Ui the conden- 
sate fraction goes to zero after the quench, whereas the 
corresponding equilibrium state is still superfluid. The 
dynamical phase diagram in Fig. [5}d summarizes our 
analysis for all kinds of quantum quenches [? ]. Let 
us finally address the changes in the dynamical behavior 
when one quenches for non integer filling. Since the 'ab- 
sorbing' Mott state disappears for n 7^ 1, it is natural to 
expect, as indeed we find, that going away from n = 1 
the dynamical transition disappears too, and transmutes 
into a cross-over that becomes more and more sharp ap- 
proaching integer filling. Overall, the resulting physical 
picture is extremely similar to the one obtained recently 
for the fermionic Hubbard model by a time dependent 
Gutzwiller approximation [10|. 

Clearly, a natural question is how much these results 
depend on the constraint of maximum two bosons per 
site. A complete analysis with an arbitrary number of 
bosons TT-b is very involved. The mapping to a classical 
system works also in these case. The classical degrees of 
freedom are the first ni, — l fractions xq, xi, . . . , Xn^-i and 



their associated canonical momenta. Unlike in the case 
rib — 2 where the classical motion is one dimensional, 
these trajectories are no longer necessarily periodic. In 
order to study their regularity we have computed numer- 
ically for rif, = 3 the largest Lyapunov exponent A [2J| of 
several trajectories. In this case a;i,a;2 are the classical 
variables and the expression of the Hamiltonian can be 
found in [l8| . Depending on the initial condition we find 
large values (A > 0.1) characteristic of chaotic trajecto- 
ries for large quenches, and small, possibly zero, values 
characteristic of periodic or quasi periodic trajectories for 
small quenches. We find again a dynamical transition, 
for n = 1 and n = 2, which are the only filling for which 
the Mott ground state exists. At the transition line, the 
trajectories are chaotic. As in the previous case, the dy- 
namical transition corresponds to a change in the form 
of the phase space trajectories: for Uf > Uf the momen- 
tum 2pi — p2 becomes unbounded, see Fig. [4^. The time 
evolution of the Xi (t) is also similar and characterized by 
oscillations that take place on longer timescales close to 
the transition. Moreover, the qualitative evolution of the 
time averaged I'I'oP (and also Xi) with Uf for a given Ui 
resembles very much the one for rif, = 2. We have also 
analyzed higher values of ri), up to nj, = 5 finding quali- 
tatively and quantitatively similar results. Actually, the 
evolution of (|\l/oP) depends very little on Ub for Ub > 2 as 
shown in Fig. |4)3 (the two curves n;, = 4 and rih = 5 differ 
by less than 0.01%). The only issue that remains open 
is the form of the singularity at the dynamical transition 
for Ub > 2. Numerical solutions of the Newton equation 
are not precise enough to answer this question. Even in 
the case of two bosons per site, for which we know that 
I'I'oP = at the transition and the singularity is log- 
arithmic, numerics alone would not be conclusive. For 
nb = 2 the singularity was due to the fact that trajecto- 
ries spend most of the time close to the Mott state. For 
Ub > 2, it is not clear whether trajectories starting with 



the same energy as the Mott state (energy zero) have to 
go arbitrary close to it. Assuming that the classical dy- 
namics is completely ergodic on the E = hypersurface, 
one could argue that this should be the case. However, 
even in this case, time averages would not coincide with 
averages in the Mott state unless the recurrence time is of 
the same order as the trapping time, a difficult question 
to address. The conclusion of the analysis performed for 
higher number of bosons is that the results for nf, = 2 
are robust and expected to hold also for the BHM with 
an arbitrary number of bosons per site, except possibly 
the form of the singularity of |\l/oP (and of the other ob- 
servables) . 

Let us now discuss the implications and the limitations 
of our findings. Clearly, the infinite dimensional limit 
neglects important dynamical and spatial fluctuations. 
This is manifest from the non damped oscillatory evolu- 
tion in time of the observables and the absence of ther- 
malization. Certainly 1/d corrections must be taken into 
account to lead to decoherence and thermalization. Nev- 
ertheless, we expect that our mean field approach should 
be able to qualitatively account for the short time dynam- 
ical behavior. As a consequence, the dynamical transi- 
tion we find should transmute into a cross-over in the 
short time dynamics for finite dimensional systems. In- 
deed results obtained for the one dimensional BHM seem 
to be in agreement with our findings [2l|. Moreover, 
we expect the dynamical transition we found to be quite 
general, at least within mean field treatments of the off- 
equilibrium dynamics. Actually, it is qualitatively iden- 
tical to the one found for the fermionic Hubbard model 
within the Gutzwiller approximation [lO| and very simi- 
lar to the one obtained by out of equilibrium dynamical 
mean field theory [9[ , where some dynamical fluctuations 
are taken into account. 

Including spatial and dynamical fluctuations would allow 
one to go beyond our mean field treatment. A good de- 
scription of decoherence and thermalization for the BHM 
could be obtained in the future within a real time gener- 
alization of the equilibrium bosonic DMFT [H, [H] . We 
would like to thank C. Kollath and M. Schiro for useful 
discussions. GB acknowledges partial financial support 
from ANR FAMOUS. 
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